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We present a detailed analysis of a model for the synchronization of nonlinear oscillators due to 
reactive coupling and nonlinear frequency pulling. We study the model for the mean field case of all- 
to-all coupling, deriving results for the initial onset of synchronization as the coupling or nonlinearity 
increase, and conditions for the existence of the completely synchronized state when all the oscillators 
evolve with the same frequency. Explicit results are derived for Lorentzian, triangular, and top-hat 
distributions of oscillator frequencies. Numerical simulations are used to construct complete phase 
diagrams for these distributions. 
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I. INTRODUCTION 

Increasingly the collective behavior displayed by 
groups of interacting dynamical components, each of 
which may be capable of a full range of complex dynam- 
ics, is essential to understanding and ultimately design- 
ing systems. Examples in biology, physics, and engineer- 
ing are diverse, ranging from understanding sensory per- 
ception to the design of antennas capable of simultane- 
ously sending and receiving signals at the same frequency. 
While in general the dynamical components may be func- 
tionally distinct and heterogeneous, in many examples 
one is interested in the collective behavior displayed by a 
group of similar coupled elements. One commonly stud- 
ied example is the synchronization of oscillating subsys- 
tems that interact. 

The ability of a group of coupled oscillators to syn- 
chronize despite a distribution in individual frequencies 
is a broadly applicable phenomenon. In physical sys- 
tems coherent oscillations can be used to enhance detec- 
tor sensitivity or increase the intensity of a power source. 
Synchronization is also important in biological systems 
where the collective behaviors in populations of animals, 
such as the synchronized flashing of fireflies or the coher- 
ent oscillations observed in the brain, serve as examples. 
Including a distribution in individual frequencies in an 
otherwise homogeneous ensemble captures some of the 
inevitable differences that group members will possess. 

Although synchronization is often put forward as an 
example of the importance of understanding nonlinear 
phenomena, the intuition for it, and indeed the subse- 
quent mathematical discussion, often reduces to simple 



linear ideas. For example, the famous example of Huy- 
gens' clocks [1] can be understood in terms of a linear cou- 
pling of the two pendulums through the common mount- 
ing support. Then the larger damping of the symmet- 
ric mode (coming from the larger, dissipative motion of 
the common support) compared with the antisymmet- 
ric mode leads at long times to a synchronized state of 
the two pendulums oscillating in antiphase. The non- 
linearity in the system is simply present in the individ- 
ual motion of each pendulum; specifically in the mecha- 
nism to sustain the oscillations. Without the drive, the 
oscillators would still become synchronized through the 
faster decay of the even mode, albeit in a slowly decay- 
ing state. Rather than this mode-dependent dissipation 
mechanism, one might expect synchronization to arise 
from the intrinsically nonlinear effect of the frequency 
pulling of one oscillator by another. Furthermore, the 
model describing the two pendulums, as well as most 
other models used to show synchronization, has dissipa- 
tive coupling between the oscillators. In contrast, many 
physical situations have mainly reactive coupling. 



Nanoscale mechanical oscillator arrays are one example 
that offer significant potential in a range of technologies 
[2], [3]. Due to the scales and numbers in these arrays, 
active control of individual oscillators poses a number 
of challenging problems. Synchronization offers a poten- 
tially appealing alternative in some applications. One no- 
table characteristic of these arrays is the predominantly 
reactive coupling due to elastic or electrostatic interac- 
tions. This context is the motivation for our work. In 
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particular, we study the model 

%m — ^(tU m Qz\z m \ )z m -(-(l \z m \ )z m H ^ ^ ^ (^n ^m) 

71=1 

(1) 

where z m is a complex number representing the ampli- 
tude r m and phase 9 m of the mth oscillator z m = r m e m . 
The natural frequency ui m of each oscillator is chosen 
from a distribution g(o->). We call the width of the dis- 
tribution w. Coefficients of the terms for nonlinear fre- 
quency pulling (a) , dissipative coupling (K ) , and reactive 
coupling (/3) serve as model parameters. The coupling is 
taken to be infinite range or all-to-all. 

In this paper we study the roles of nonlinear frequency 
pulling and reactive coupling in the absence of the dissi- 
pative coupling using Eq. (1) with a ^ 0, (3 ^ 0, K = 0. 
In addition the distribution g(uj) of the w m must be spec- 
ified. We study the case of positive a and /3; for a sym- 
metric distribution g(tu) the results are the same chang- 
ing the sign of both a and (3. 

Alternatively, when only nonlinear saturation and dis- 
sipative coupling are present (a = f3 = 0, K ^ 0) Eq. (1) 
reduces to 

K N 

Z m = (iw m + 1 - \z m \ 2 )z m + J^^i^n - Zm)- (2) 

n— 1 

The behavior for general w and K of Eq. (2) has been 
analyzed by Matthews et al. [4]. If the width w of the 
distribution g(uj) is narrow, so that the time evolution 
of the magnitudes \z m \ is fast compared with that of the 
phase dispersion, and K is small, \z m \ rapidly relaxes to 
a value close to unity, and the only remaining variable 
for each array member is the phase 8 m . Equation (2) 
can then be reduced to a simple form [5] , [6] , often known 
as the Kuramoto model 

is 

6 m = ui m + — 2J sin(6»„ - 6 m ), (3) 

71=1 

that has been the subject of numerous studies [7]. In the 
absence of coupling each oscillator in this model would 
simply advance at a rate that is constant in time, but 
with some dispersion of frequencies over the different ele- 
ments. Equation (3) is an abstraction from the equations 
describing most real oscillator systems, leaving out many 
important physical features. 

The complex notation of Eq. (1) suggests the intro- 
duction of a complex order parameter tp to measure the 
coherence of the oscillations 

1 N 

y J = Re^ = -J2r n e^, (4) 

71=1 

with r„ = 1 for the Kuramoto model, and then the equa- 
tion of motion can be written [take the imaginary part 
of tpe' l6m in Eq. (4) to evaluate the sum appearing in 
Eq. (3)] as 

m =w m + KRsm(Q-6 m ). (5) 



Thus the behavior of each oscillator is given by its ten- 
dency to lock to the phase of the order parameter. The 
term KRsin(Q — 6 m ) acts as a locking force, and lock- 
' ing occurs for all oscillators with frequencies satisfying 
\uj m \ < KR, with the locked oscillator phase given by 
O + sin (u> m /KR). The magnitude R of the order pa- 
rameter must then be determined self-consistently via 
Eq. (4). The generalization of the locking force to ap- 
ply to the model Eq. (1) will be a conspicuous feature of 
our work. 

Equation (3) is known to show rich behavior, including, 
in the large N limit, a sharp synchronization transition 
at some value of the coupling constant K = K c [6] , which 
depends on the frequency distribution of the uncoupled 
oscillators g(w). The transition is from an unsynchro- 
nized state with if) = in which the oscillators run at 
their individual frequencies, to a synchronized state with 
tjj 7^ in which a finite fraction of the oscillators lock to 
a single frequency. The transition at K c has many of the 
features of a second order phase transition, with univer- 
sal power laws and critical slowing down [6], as well as a 
diverging response to an applied force [8]. 

We present a detailed analysis of the model Eq. (1) 
with K = describing the synchronization of an oscilla- 
tor ensemble involving reactive coupling between the ele- 
ments, which then leads to synchronization through non- 
linear frequency pulling. We begin in Sec. II by deriving 
Eq. (1) as a description of arrays of nanoelectromechan- 
ical oscillators. We then discuss the solutions to Eq. (1) 
for a variety of symmetric distributions in intrinsic fre- 
quencies g(iu). Common solution types exist for the three 
distinct types of frequency distributions we studied. In 
Sec. Ill we introduce these solutions and the measures we 
use to describe them. We begin the analysis in Sec. IV 
by moving to a continuum description to derive synchro- 
nization conditions. We are able to analyze the existence 
of two behavior types in closed form: the unsynchro- 
nized solution in Sec. V and the fully locked synchro- 
nized behavior in Sec. VI. Results from our analytical 
arguments are combined with those from simulations to 
present the solutions and associated phase diagram for 
each frequency distribution in Sec. VII. We begin these 
by considering an unbounded symmetric frequency dis- 
tribution in the form of a Lorentzian in Sec. VII A. Con- 
tinuing in this section we then approximately double the 
distribution width to discuss some interesting character- 
istics of the well-ordered synchronized solutions as well 
as changes to the bifurcations. In contrast, Sec. VII B 
presents the results for a bounded distribution when all 
frequencies over some range are equally likely (a top-hat 
distribution). This is followed in Sec. VII C by the case 
of a symmetric unimodal frequency spread, namely a tri- 
angular distribution. Finally, conclusions are made in 
Sec. VIII. 

A brief account of some results for the model Eq. (1) 
has been reported previously [9]. 
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II. CONNECTION WITH MECHANICAL 
OSCILLATORS 

Although the main focus of this paper is analyzing the 
behavior of Eq. (1) we first show how such an equation 
might arise from the equations of motion of physical os- 
cillators. As an example consider the system defined by 

x n +^x n -v{l-x 2 n )x n +ax\-D[x n ~\{x n+1 +x n -i)] = 0. 

(6) 

The first two terms describe uncoupled harmonic oscilla- 
tors. We suppose the uncoupled oscillators have a linear 
frequency that is near unity by an appropriate scaling of 
the time variable 

ujI = 1 + S„ , 5 n < 1. (7) 

The third term is a negative linear damping, which rep- 
resents an energy source to sustain the oscillations, and 
positive nonlinear damping, so that the oscillation am- 
plitude saturates at a finite value. Again this satura- 
tion value is chosen to be of order unity by an appropri- 
ate scaling of the displacements x n . For an example of 
an effective negative linear damping term in a microme- 
chanical oscillator see references [10, 11]. One could also 
imagine implementing such an effect with an electronic 
feedback loop sensing each oscillator velocity and driv- 
ing the oscillator with an appropriate phase. The first 
three terms of Eq. (6) comprise a set of uncoupled van 
der Pohl oscillators. The term ax\ describes a stiffening 
of the spring constant (for a > 0) and is a reactive non- 
linear term that leads to an amplitude dependent shift 
of the resonant frequency. With v = this would give 
us Duffing oscillators. The final term is a nearest neigh- 
bor coupling between the oscillators, depending on the 
difference of the displacements. This is a reactive term, 
that conserves the energy of the system. Others [12] have 

I 



The collection of terms proportional to e lt on the right- 
hand side of Eq. (11), called the secular terms, act like a 
force driving the simple harmonic oscillator on the left- 
hand side at its resonance frequency. The sum of all the 
secular terms must vanish so that the perturbative cor- 
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considered nonlinear oscillators coupled through their ve- 
locities; this is a dissipative coupling that would lead to 
K 7^ in the amplitude-phase reduction. 

The complex amplitude equation (1) holds if the pa- 
rameters a, -D, 8 n are sufficiently small. In this case the 
equations of motion are dominated by the terms describ- 
ing simple harmonic oscillators at frequency one, and the 
time dependence remains close to e ±rf . To formalize the 
smallness of u,a,D,8 n we introduce a small parameter 
e and write v = sD, S n = sS n , a = ea, D = eD. The 
oscillating displacement is written as a slow modulation 
of oscillations at frequency one, plus corrections [13] 

x n (t) = [A n (T)e H + c.c] + ex\l\t) + ... (8) 

with T = et a slow time scale [14]. The variation of 
A n (T) gives us the extra freedom to eliminate secular 
terms and ensure that the perturbative correction xffl (i), 
as well as all higher-order corrections to the linear re- 
sponse, do not diverge (as they do if one uses naive per- 
turbation theory). Using the relation 

• dA n dA n , 

. An =^r = e iT= eA ^ . (9) 

(denoting a time derivative with respect to the slow time 
T by a prime) we calculate the time derivatives of the 
trial solution (8) 

x n = ([iA n + eA' n ]e lt + c.c.) + ex^(t) + ... (10a) 

x n = ([-A n + 2ieA' n + e 2 <]e lt + c.c.) + + . . . 

(10b) 

Substituting these expressions back into the equation of 
motion (6), and picking out all terms of order e, we get 
the following equation for the first perturbative correc- 
tion 



(11) 

I 

rection Xn (t) will not diverge. (Terms varying as e ±3lt 

contribute a finite response to x n .) This gives us an 
equation for determining the slowly varying amplitudes 
A n (T) 



xlP + 4 1J = -SnA n - (2iA' n e lt + c.c.) + v (iA n e lt + c.c.) (1 - (A n e u + c.c.) 2 ) 
- a (A n e u + c.c.) 3 - D[(A n - \{A n + A n ^)) e lt + c.c.}. 



2A' n = {v + iS n )A n -(v- Ma) \A n \ 2 A n - iD[A n - \{A n+1 + A„_i)] = 0. (12) 
More informally, we might write x n (t) = z n (t)e lt + c.c + ■ ■ ■ so that A n (T) — > z n (t) and Eq. (12) can be written in 
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terms of the original "small" parameters without the formal scaling 

2i„ = (v + iS n )z n - (v - 3ia) \z n \ 2 z n - iD[z n - \{z n +i + z n -i)) = 

I 



(13) 



With a rescaling of time t = i/t/2 Eq. (13) reduces to the 
form Eq. (1) except that in Eq. (1) the nearest neighbor 
coupling is replaced by the all-to-all coupling convenient 
for theoretical analysis. 



III. EXAMPLES OF DYNAMICAL STATES 

In this section we introduce the types of dynamical 
states encountered for the model Eq. (1) as well as di- 
agnostic tools to characterize these states. The types 
of states we find for the different distributions investi- 
gated are largely the same, and so we use a particular 
example — a Lorentzian distribution of oscillator frequen- 
cies with some convenient choice of parameters a and 
(3 — and discuss how the behavior depends systematically 
on the distribution and other parameters later in the pa- 
per. We are mainly interested in the behavior for large 
numbers of oscillators TV — > oo. For the numerical simu- 
lations we are of course restricted to finite TV (we typically 
use TV = 1000, but have gone up to TV = 100, 000 to in- 
vestigate some subtle effects) . In our discussion we focus 
on those results that we expect to be largely independent 
of TV for large TV. 

A key diagnostic for synchronization is the complex 
order parameter \I/(i) defined by Eq. (4), introduced by 
Kuramoto [6], with magnitude R and phase where 
r n e lSn = z n . In the large TV limit, we could use a nonzero 
value of the order parameter at each time as the basic cri- 
terion for a synchronized state. A synchronized periodic 
state with frequency f2 would then have = Re iQt 
with R, Q constants. More complicated dynamical states 
might also occur. For a finite number of oscillators TV a 
precise diagnostic is harder, since there are fluctuations of 
order TV -1 / 2 that make the instantaneous \& nonzero even 
in an unsynchronized state. For a synchronized periodic 
state we could require that the time average ^e _lf2 ( JV )*\I/) 
scales as TV as the number of oscillators changes for some 
frequency fi(TV) that becomes constant for large enough 
TV, and use this as a measure of the synchronization. In 
practice we use the simpler criterion that the time aver- 
aged magnitude (R(t)) t = (\ip{t)\) t 1S nonzero and does 
not appear to decrease to zero as TV increases. (A ^ fluc- 
tuating about zero with an amplitude of order TV -1 / 2 will 
of course lead to a nonzero (R(t)) t of this same order.) 
This definition can also be applied to aperiodic synchro- 
nized states. 

Another useful diagnostic uses the actual frequency of 
each oscillator defined by 



where t is some long averaging time and to is a start- 
ing time sufficient to eliminate transients. A frequency 
locked state has a fraction of oscillators with the same 
frequency u> n . (The fraction should be 0(1) and a value 
not decreasing to zero as TV increases). If not all the os- 
cillators have the same frequency (i.e. the fraction is not 
unity) we call the state partially locked. If all the oscil- 
lators have the same average frequency, we call the state 
fully locked. To make contact with the analytic results 
we actually use a stricter criterion, and also require the 
magnitudes r n to be time independent in the fully locked 
state. Usually we find that a nonzero (R(t)) t is associ- 
ated with frequency locking, but this is not always the 
case. 

Phase locking is a stricter requirement than Eq. (14). 
For phase locking we would require that 8 n (t) — m {t) 
does not diverge as t — ■+ oo for to, n taken from some 
finite fraction of the oscillators. (Frequency locking is, 
for example, consistent with phase differences that grow 
diffusively proportional to i 1 / 2 .) We do not investigate 
this stricter locking criterion. 

To investigate whether each locked oscillator in a fre- 
quency locked state is tightly locked in phase to the phase 
of the order parameter, or is fluctuating about this value, 
we define a phase synchronization index for each oscilla- 
tor 



Xn — 1 „ 



iO n (t) 



(15) 



with 9 n the phase of the nth oscillator relative to the 
order parameter phase 



e n (t) = e n (t) - e(t). 



(16) 



When the phase of oscillator n is locked to the order 
parameter one, 9 n {t) is constant and the phase synchro- 
nization index is unity, whereas as 6 n (t) tends towards a 
uniform distribution from to 27r the index approaches 
zero. We define the number of oscillators with \n very 
close to one as the tightly locked cluster. 

We now show results of numerical simulations for the 
Lorentzian distribution of frequencies. We have chosen to 
cut off the distribution at some large oj c otherwise there 
would be a few very fast oscillators restricting the time 
stepping of the numerics. Thus we use 



9(0) 




for I a; I < uj c 
for \u>\ > oj c 



(17) 



for a given choice of g(0), with w then fixed by the nor- 
malization condition 



■jj n 



[o n (to + t) -e n {t )} 



(14) 



2iu<7(0)tan (u) c /w) = 1. 



(18) 
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FIG. 1: Distribution g(u>) of frequencies for cutoff Lorentzian 
with g(0) = 1 and cutoff w c = 8 for N = 1000 oscillators. The 
main graph shows the individual frequencies, and the inset a 
histogram. 



(It is useful to parameterize the distribution in terms 
of g(0), since this quantity determines the values of a 
and /3 at which synchronization occurs in the large a(3 
limit where the model reduces to the Kuramoto phase 
model Eq. (3).) In presenting the results we choose a 
distribution width such that g(0) = 1. For a cutoff lo c = 8 
this gives w ~ (/tt)~\ with / = 0.974. With no cutoff, 
we would have / = 1. The distribution of frequencies is 
generated from a uniform distribution of N values v n on 
the interval — \ < v„ < | using the transformation 

uj n = u;tan(7r/u„). (19) 

We present results for a deterministic set of v n 



n-\(N + l) 



/(N-l), 



(20) 



but have also used v n generated randomly on the interval. 
The distribution of frequencies for N = 1000 oscillators 
and a cutoff lu c = 8 is shown in Fig. 1. 

A plot of the dependence of the mean magnitude of 
the order parameter (R) t as a function of f3 for fixed 
value of a = 1.5 and N = 1000 is shown in Fig. 2. This 
plot shows three types of states: an unsynchronizcd state 
with (R) t essentially zero; and two synchronized states, 
one with small (R) t , growing from to about 0.3 as (3 
decreases below about 5.2, and one with large (R) t that 
exists for all /3 > 0.6. (We will discuss later such issues 
as whether the growth of the large (R) t is continuous or 
discontinuous near /3 = 0.6.) 

The variation of R(t) for representative examples of 
these three states is shown in Fig. 3. The lower trace 
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FIG. 2: Solutions observed in simulations of TV = 1000 os- 
cillators with the cutoff Lorentzian distribution of intrinsic 
frequencies of Fig. 1. Shown is the time averaged order pa- 
rameter magnitude {R)t over a range of values at constant 
a — 1.5. Solid lines are the observed solutions. The over- 
lapping symbols are representative results from simulations 
following different solution branches with increasing j3 (A) 
and decreasing j3 (□ and O)- Time traces associated with 
the solid symbols are shown in Fig. 3. The dot-dashed line 
is at a value of {R)t — N^ 1 ^ 2 — 0.032, the order of magni- 
tude expected for random fluctuations of the order parameter 
about zero. 



shows the unsynchronized state at (3 = 6.0 in Fig. 2. The 
fluctuations in R(t) are consistent with fluctuations of 
the order parameter about zero with magnitude of order 
N^ 1 / 2 as expected for a collection of N oscillators with 
different frequencies and random phases. The average 
frequency distribution uj n (not shown) is unchanged from 
the bare distribution u> n shifted by a + (3. This shift can 
be understood as arising from the nonlinear frequency 
shift with \z n \ = 1, and the coupling to a distribution of 
oscillators with \z m \ = 1 and random phases. 

The state corresponding to the middle trace in Fig. 3 
is the low amplitude state at j3 = 4.0 in Fig. 2. The 
mean order parameter magnitude (R) t — 0.17 is much 
larger than N~ x / 2 ~ .032, suggesting this is a synchro- 
nized state. The corresponding frequency distribution 
and synchronization index are shown in Fig. 4. The fre- 
quency distribution Q n (a) shows a small plateau of con- 
stant frequency over about 60 oscillators towards the high 
frequency end of the distribution; these are the frequency 
locked oscillators. The synchronization index (b) shows 
that Xn approaches unity for most of the locked oscilla- 
tors, which means that z n for these oscillators is essen- 
tially time independent once the rotation of the phase at 
the order parameter frequency is subtracted out. A care- 
ful scrutiny of the two panels reveals that the plateau 
in Co n is sharper and more extended than the one in Xm 
so that not all the frequency locked oscillators are time 
independent in the rotating frame. The dynamical state 
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FIG. 3: Order parameter magnitude 7? as a function of ti 
observed at the solid symbols in Fig. 2. The three traces 
for the unsynchronized state at /3 = 6.0 (lower trace), and the 
small amplitude synchronized state (middle trace) and large 
amplitude synchronized state (upper trace) both at /3 = 4.0. 




oscillator 

FIG. 4: Frequency distribution and synchronization index for 
the small amplitude synchronized state in Fig. 3: (a) - actual 
frequency distribution (solid line) and bare frequency distri- 
bution shifted by a + j3 (dotted line); (b) - synchronization 
index Xn- 



is actually quite complicated, as a review of z n (t) shows. 

Figure 5 shows a plot of a snapshot of the complex 
amplitude z n of each oscillator. The "+" is the order 
parameter. The other points correspond to z n for each 
oscillator. It is useful to study the dynamics of this plot 
after rotation at the mean order parameter frequency is 
eliminated. As time evolves, the square symbols remain 
fixed in such a plot (except for very small fluctuations): 
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FIG. 5: Snapshot of z„ and order parameter <i? in the com- 
plex plane at a time after transients have decayed. The "+" 
towards the center of the plot gives the value of the order 
parameter. The squares are oscillators locked to the order 
parameter frequency (i.e. those with u) n = Q): solid squares 
are stationary when the e' nt dependence is removed, open 
squares show an additional dynamics rotating around the tail 
of the fixed distribution in the rotating system. The "x" are 
oscillators that are not locked to the order parameter. 



these represent oscillators that are locked to the order 
parameter. For the solid squares the complex ampli- 
tudes are essentially time independent once the phase 
of the order parameter is extracted. The oscillators rep- 
resented by an "x" on the other hand rotate clockwise or 
anticlockwise about the origin: these correspond to un- 
locked or running oscillators. The open squares exhibit a 
more complicated dynamics undergoing small amplitude 
orbits around the tail of the locked oscillator distribution. 
These oscillators arc locked to the order parameter, since 
the difference of their phases from the order parameter 
phase does not drift over arbitrary long times. The val- 
ues of uj n for these oscillators are on the locked plateau 
in Fig. 4. However the amplitudes are not constant in 
the rotating plot, the values of Xn are less than unity, 
and the oscillators contribute to fluctuations of the order 
parameter. Thus the fluctuations of R(t) shown in Fig. 
3 are not just due to finite TV effects, and we believe they 
would persist in the ./V — > oo limit. We have not explored 
larger TV to pin down whether these intrinsic fluctuations 
are periodic or aperiodic in the large TV limit. 

As we discuss in more detail later, for a bounded distri- 
bution of frequencies it is possible to find a low amplitude 
synchronized state with R 0, but with a smooth fre- 
quency distribution showing that there is no frequency 
locking. In fact for this state the distribution of actual 
oscillator frequencies does not overlap the order parame- 
ter frequency. The nonzero order parameter is caused by 
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a systematic slowing of the phase rotation of the oscilla- 
tors in the vicinity of the order parameter phase, rather 
than by a fraction of the phases becoming locked to the 
order parameter phase. For an unbounded distribution 
such as the Lorentzian, this same effect occurs but is 
supplemented by the small fraction of locked oscillators. 

10 1 1 1 1 1 1 1 1 1 1 1 
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FIG. 6: Same as in Fig. 4 but for the large amplitude syn- 
chronized state. 



For the same parameter values as for the previ- 
ously described low amplitude synchronized state for the 
Lorentzian distribution there is a second dynamical state 
that may be reached depending on the initial conditions 
of the simulation. This state has a much larger value of 
the order parameter (R) t — 0.9 and the plot of the dis- 
tribution of frequencies, Fig. 6(a), shows a correspond- 
ingly larger fraction of locked oscillators. (This should 
be compared with panel (a) of Fig. 4.) For the cutoff 
Lorentzian distribution used, all the oscillators at the 
high frequency end of the distribution arc locked, leaving 
only a small fraction of free running oscillators at the low 
end of the frequency distribution. For the N — ► oo limit 
and a Lorentzian distribution without a cutoff, unlocked 
oscillators remain at both the high and low frequency 
ends of the distribution. For a bounded distribution of 
bare frequencies, such as triangular or top-hat, a fully 
locked state in which all the oscillators rotate with the 
same frequency and fixed magnitude may be found. 

IV. FORMULATION OF SYNCHRONIZATION 

We now turn to the analysis of the synchronization of 
oscillators described by Eq. (1). Since we are interested 
in the behavior for a large number of oscillators, it is 
convenient to go to a continuum description, where we 



label the oscillators by their uncoupled linear frequency 
u> = ui n rather than the index n, z„ — ► z(u>). Introducing 
the order parameter Eq. (4), the oscillator equations can 
be written in magnitude-phase form as 

d t 6 = Lj + a(l-r 2 ) + — cosfl (21a) 
r 

d t r = (l-r 2 )r + /3Rsin6 (21b) 

where 8 = 8 — O is the oscillator phase relative to that of 
the order parameter as before, and ui is the bare oscillator 
frequency shifted by a + f3 and measured relative to the 
order parameter frequency fi = 

w=w-a-P-Cl. (22) 

Note that if the order parameter is zero R = 0, the mag- 
nitude r will relax to 1 , and the nth oscillator will evolve 
at the frequency u>„ — a — f3. The first component of this 
frequency shift is just the nonlinear shift at r = 1, and 
the second is from the interaction of each oscillator with 
the incoherent motion of the other oscillators. 

At each time t the oscillators are specified by the distri- 
bution p(r, 8, u>, t) where p{r, 8, cD, t)r dr dd is the fraction 
of the oscillators with shifted frequency Co that at time t 
have magnitude between r and r + dr and shifted phase 
between 8 and 8 + dd. The order parameter is given by 
the self-consistency condition 

R = (re i3 ^ = J dCog{Lj) J r dr d8p(r,8,Co,t)re lS . (23) 

where g(u>) is the distribution of oscillator frequencies 
expressed in terms of the shifted frequency u>. It is useful 
to split this expression into real and imaginary parts. 
The imaginary part is 

J duigluj) J r dr d6 p(r, 9, u),t)r sin 6 = 0. (24) 

Because the phases 8 are measured relative to the orien- 
tation of the order parameter, this expresses the fact that 
the components of the complex amplitudes r sin 8 normal 
to the order parameter must average to zero. Note that 
unlike the cases of the Kuramoto model [Eq. (3)] and the 
dissipativcly coupled complex amplitude model [Eq. (2)] 
studied by Matthews et al. , this condition is not trivially 
satisfied even for the case of a symmetric distribution 
g{ijj), and in fact serves to determine the frequency fl of 
the order parameter. The frequency f2 is also not trivially 
related to the mean frequency of the oscillator distribu- 
tion. The real part of Eq. (23) is 

J duig{ui) J r dr d8p(r, 8, w, t)r cos 6 = R. (25) 

This is the condition that the components r cos 8 along 
the direction of the order parameter must average to the 
magnitude R. This condition serves to sclf-consistently 
fix the value of R. 
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The expectation that Eqs. (21) might lead to synchro- 
nization follows from the behavior for narrow frequency 
distributions and large a. If the width of the distribu- 
tion of frequencies is small compared to the relaxation 
rate of the magnitude, which is of order one in the time 
units used in Eq. (1), the magnitude relaxes rapidly to 
the value given by the instantaneous value of 9, i.e. to 
the solution of 



(1 



= -/3iisin( 



(26) 



If r is close to one, which we will sec applies at the onset 
of synchronization for large a, this gives 



1 



-/3RsmO. 



(27) 



In this case Eq. (21a) becomes, ignoring (3R compared 
with a[3R, 



dt9 ~ to — a(3R sin 9. 



(28) 



This equation is the same as the one derived from the 
Kuramoto model in Eq. (5) with a(3 playing the role of 
the coupling constant K, and therefore predicts an onset 
of synchronization at a(3 = 2(irg(0))~ 1 [6]. 

To uncover more fully the behavior of Eq. (1) we con- 
sider two issues: the onset of synchronization, detected as 
the linear instability of the unsynchronized R = state, 
as a function of a and P for some given frequency distri- 
bution g(ui); and the existence of a fully locked state for 
large values of a/3. 



V. ONSET OF SYNCHRONIZATION 

We first consider the initial onset of partial synchro- 
nization from the unsynchronized state in which each os- 
cillator runs at its own frequency w„, We identify syn- 
chronization through a nonzero value of the order pa- 
rameter. This often arises because a finite fraction of the 
oscillators become locked to the same frequency, which 
we call a partially (fully if the fraction is unity) locked 
state. However we also find situations where R ^ 0, but 
there is no frequency locking. 

The onset of synchronization can be determined by 
a linear instability analysis of the unsynchronized state. 
This is calculated by linearizing the distribution p around 
the unsynchronized distribution, which is a uniform 
phase distribution at r = 1, and seeking the parameter 
values at which deviations from the uniform phase dis- 
tribution begin to grow exponentially. This follows the 
method of Matthews et al. [4] , although care is needed 
in the analysis due to the more important role the mag- 
nitude perturbations r play in the present case. 

Introducing the small expansion parameter e charac- 
terizing the small deviations from the unsynchronized 
state, we write 



p(r,6,Q,t) = {2nr)- l 8[r-l- 



■erx(8,u),t)}[l+eh{6,u>,t)\. 

(29) 



Note that for e = this does indeed give the appro- 
priately normalized distribution for the unsynchronized 
state, in which all the oscillators have unit magnitude, 
r = 1, and the phase distribution is constant. Also, for 
nonzero e, p remains normalized to linear order 



dr 



d9rp(r,9,Cd,t) = l + 0(e 2 ), 



(30) 



providing the average of f\ over 9 is zero. 

The equation for the evolution of the radial perturba- 
tion of the distribution n is given by noting that 



dr dr - dr 

~dt = W t+ m- 



(31) 



The left hand side is evaluated by expanding Eq. (21b) 
to first order in e and also expanding the magnitude of 
the order parameter R = eR\ + • • • . , and the right hand 
side by the replacement d t 9 = Co + 0(e) and assuming an 
exponential growth or decay of the perturbation dr/dt = 
edri jdt = \er\ . The result is 



dri 
~d~9 



u + (X + 2)n = PRisint 



Equation (32) is solved by 

n = i?i(Acos6» + Bsm9), 

with 

A- a g 

P Co 2 + (X + 2) 2 ' 

R _o (A + 2) 
P Cd 2 + (X + 2f 



(32) 
(33) 

(34a) 
(34b) 



To extract the equation for /i integrate the equation 
for the conservation of probability 



dp 

dt 



V • (pv) = 



(35) 



over radius. Here v is the velocity in complex amplitude 
space, which in polar coordinates is (dtr, rdt9) given by 
Eqs. (21). Again replacing dfi/dt by A/i, and evaluating 
ri from Eqs. (33,34), this gives at 0(e) 

Xfi+udgfr = 2ai?i(-A sin 6+B cos fl)+#Ri sin£. (36) 

This is solved by 

h = Ri(C cos 9 + D sin0) (37) 

with 

2a(A 2 + 2A - Co 2 ) - up 2 + (A + 2) 2 ] 



C = p 
D = 



(Co 2 + A 2 ) [Co 2 + (A + 2) 2 ] 
4a^(A + 1) + \[Co 2 + (A + 2) 2 ] 
(Co 2 + A 2 ) [Co 2 + (A + 2) 2 ] 



(38a) 
(38b) 
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We now evaluate the self-consistency condition 
Eq. (23) to first order in e. The imaginary part is 

J dujg{u>) J d9(l + er 1 + ---)(l + ef 1 + ---)sm9 = 0, 

(39) 

which to first order in e gives 

J dCog(Co)[B + D] =0. (40) 
Similarly the real part of the self-consistency condition is 

J du> g{w) [A + C] = 2. (41) 

We want to evaluate Eqs. (40) and (41) at the onset 
of instability where the growth rate A — > 0. We can set 
A = in Eqs. (40,41) with Eqs. (34) and (38) except 
in terms with Co 2 + A 2 in the denominator, since such 
terms may give large contributions to the integral from 
the region of small Co. A term involving just A/ {Co 2 + A 2 ) 
gives a finite integral, but if this is multiplied by powers 
of Co or A the integral goes to zero in the A — » limit. 
Similarly for a term involving Co /{Co 2 + A 2 ) we must take 
the limit A — > after doing the integral (this is equivalent 
to the principal value integral), whereas if this term is 
multiplied by powers of Co we can put A = immediately. 

The needed integrals are 



h 




^ 2 +4' 




(42a) 


h 


= lim 

A^O 


S 9{U) C0 2 + 


A 2 ' 


(42b) 


h 


"A 


^Co 2 +4> 




(42c) 


h 


= lim 

A^O 


S 9{U) C0 2 + 


A 2 = 7T.9(^ = 0). 


(42d) 



The imaginary part of the self consistency condition 
Eq. (40) becomes 

21 1 + ah - ah + h = 0, (43) 

and the real part reduces to the condition for (3 C 

fa = 2{-I 3 + ah - 2ah - h)' 1 - (44) 

We have explicitly evaluated the integrals for top-hat, 
triangular, and Lorentzian distributions of bare frequen- 
cies. These results will be presented after we discuss full 
locking. 

VI. FULL LOCKING 

We define the fully locked state as one in which all 
the phases are rotating at the same frequency as the or- 
der parameter, and the magnitudes are constant in time. 



These solutions are defined by Eq. (21a) with dtO = 0, 
which with Eq. (26) can be written 

Co = — (asin0-cos0) = F(6), (45) 

where the solution to the cubic equation (26) for r{9) is 
to be used to form the function of phase alone F{9). The 
function F{8) acts as the force pinning the locked os- 
cillators to the order parameter, generalizing the notion 
introduced below Eq. (5), and plays a central role in our 
discussion of locking. A particular oscillator, identified 
by its shifted frequency Co, will be locked to the order 
parameter if Eq. (45) has a solution 9 = F^ 1 ^) [and 
then r is given by solving Eq. (26)] and if this solution 
is stable. The stability is tested by linearizing Eqs. (21) 
about the solution. The fully locked solution will only 
exist if stable, locked solutions to Eq. (45) exist for all 
the oscillators in the distribution. We are thus led to in- 
vestigate the properties of the function F{9). In addition 
the self consistency condition Eq. (23) must be satisfied. 

For small f3R, the magnitude r{9) given by Eq. (26) re- 
mains bounded away from zero for all 8, and the function 
F{9) varies continuously between minimum and maxi- 
mum values .Fmin < F < -Finax- In this case, we im- 
mediately see that the fully locked solution only oc- 
curs for bounded distributions, g{Co) nonzero only be- 
tween finite (D m i n and <D max - In such cases we define 
w ro ax — Wmiii, which is equal to the range of unshiftcd 
frequencies o-> max — w m i n , as w the width of the distri- 
bution. More generally, although F{9) can vary over an 
infinite range (because r may become zero) we find that 
only a finite range yields stable solutions, so that again 
complete locking only occurs for a bounded distribution 
of oscillator frequencies. 

We first look at the fully locked solution for large values 
of aj3. In this case the phases of the locked oscillators 
cover a narrow range of angles, since the range of the 
pinning force F becomes large for large a/3. The imagi- 
nary part of the self consistency condition Eq. (24) shows 
that the range of phases must be around 9 = 0. Equa- 
tion (45) can now be approximated by expanding around 
9 = (note r ~ 1 here) and becomes for large a 

Co = Lo-a- /3 - ft ~ -/3R(1 - a6). (46) 

The imaginary part of the self-consistency condition re- 
duces to (9) = (the average is over the distribution of 
frequencies), and the real part to simply R ~ 1. Finally, 
averaging Eq. (46) over the distribution of frequencies 
fixes the order parameter frequency (the common fre- 
quency of all the oscillators) 

ft ~ (u) - a. (47) 

Thus the order parameter, and all the oscillators, evolve 
at a frequency given by the mean of the distribution g(to) 
shifted by the nonlinear effect for r = 1. 

We now investigate the limit to the fully locked regime 
as we lower a. We first summarize the argument, and 
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then present the details. The fully locked solutions are 
determined by a rather complicated set of interconnected 
equations. They can be found by the following algorithm. 
For fixed values of a, B = f3R solve Eq. (26) for real 
positive r(9) and hence calculate F(6). Using Eqs. (21) 
the stability of each solution is tested: the eigenvalues of 
this analysis are (see Eq. (60) below) 



l-2r^±\/l 



52 



2r 2 + 2r 4 - 2arB cos 6, (48) 



with cos 9 = ±\/l — r 2 (l — r 2 ) 2 / B 2 . From this we can 



identify a range F^) n < F < F^i x corresponding to the 
range of existence and stability of locked oscillators. A 
fully locked solution must then satisfy the constraints 



,-, ■> p( s ) . ,-, <- p(s) 
"'mm ^ ^min' "'max ^ r max 



(49) 



together with the condition given by the imaginary part 
of the self consistency condition 



g(io) r{w) sin[0(u;)] dto = 0, 



(50) 



where 6(u>) = F 1 (ui), and r(ui) is then the solution to 
Eq. (26). The boundary of the fully locked region occurs 

either when ui min = F^ n or when w max = infix- This 
equation can be interpreted as fixing the order parameter 
frequency fl in terms of a and B. Notice that Eq. (48) 
shows that for r < 1/V2 one eigenvalue certainly has a 
positive real part indicating instability, so that for stable 

solutions r > 1/ \[2 and F^ n is finite. Since w m in or w ma x 
is now determined, and w max — Lu m i n = w, Eq. (50) is an 
implicit equation relating the values of a, B and w at 
the locking transition. To complete the solution, the real 
part of the self consistency equation 



g(oj) r(oj) cos[8(lo)] doj = R 



(51) 



then serves to fix R at locking, from which the value of 
P = B/R at the transition to full locking can be found. 



A. Existence of individual oscillator locked solution 

We first consider the existence of a locked solution 
for an individual oscillator, i.e. a stationary solution of 
Eqs. (21). Equation (26) gives the cubic equation for r(9) 
for each B = (3R 



(1 - r 2 )r + B sin6> = 0. 



(52) 



Of course, for the physical solution r must be real and 
positive. Thus we need to analyze the properties of the 
real positive solutions to 



(l-r 2 )r + X = 0, 



(53) 



Roots 








CVJ 


\ 1 x 


2 







FIG. 7: Solutions to the cubic equation Eq. (53) as X varies. 
The dashed lines are y = ±l/\/3. 



as X varies. 

For X = the solutions to the cubic are r = ±1, 0. 
As \X\ increases, the solutions remain real until two of 
the roots collide and become complex. Since the sum of 
the roots to Eq. (53) is zero, at the collision the equation 
takes the form 



(r-a) 2 (r + 2a) = 0. 



(54) 



Matching coefficients shows that collision occurs at X 



± 



V27 



when r 



V3' 



The form of Eq. (53) is actually 



already in what is known as the "depressed form" of a 
cubic equation, for which the solution is relatively sim- 
ple. Inspecting the form of these solutions shows that 
for \X\ > -j= there is one real solution and a complex 
pair in the form 2a, —a ± ib. The product of the roots is 
determined by the constant in the cubic, giving 



2a(a 2 + b 2 )=X 



(55) 



and so for X < —^j=j= the real root is negative, and for 
X > -j= the real root is positive. These results are 
confirmed by numerical solution as shown in Fig. 7. 

Thus we find the following behavior for the real posi- 
tive solutions to Eq. (52). For B < 2/V27 there is a root 
that is positive with r > l/v3 for all 9, a root that varies 
between positive and negative values with magnitude less 
than 1 / v3, and a root that is negative for all 9. We have 
seen that the stability analysis shows that any solution 
with r < 1/ y/2 is unstable. Thus only the first root is rele- 
vant (see Fig. 8). For B > 2/ \/27 real positive roots exist 
for — 9b < 9 < tt+9b with 6b = sin -1 B ^f - In the range 
where there are two real positive roots, —9b < < and 
7r < 9 < 7r + 6b, only the larger may be in the stable 
range r > 1/ y/2, (see Fig. 9) . 

Having found the stationary solutions r{9), the condi- 
tion for stationary phase dt9 = can be written 



with 



uj = F{9) 



(3R 



F{9) = (a sin 6>- cos 6). 
r{6) 



(56) 
(57) 
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FIG. 8: Plot of r(0) for S < 2/V27. The dashed line is 
r = 1/V3 and the dotted line r = l/y/2. The condition 
r > l/\/2 is necessary (but not sufficient) for the solution to 
be stable. The lower branch always satisfies r < l/v3, and 
so only the upper branch may have stable regions. 



with 



5 = 1- 



B 2 

r 2 



2r 2 



2r 4 - 2a.Br cos 6 



(61) 



This immediately shows us r > 1/V2 is a necessary con- 
dition for stability, since Re A + > 1 — 2r 2 (with the equal- 
ity if S is negative so that X± are complex). 

Let us first seek the condition for a root A± to become 
zero, signaling a stationary bifurcation. It is convenient 
to go back to the original equations Eqs. (21) in the form 



dt§ = u-f(r, 9) 

d t r = (1 -r 2 )r + l3R sm.9 



with 



f(r,6) = -a(l-r 2 ) 



0R 



cos 9. 



(62a) 
(62b) 



(63) 



r 

y 























1 



e/ji 



FIG. 9: Plot of r(6) for B > 2/V27. Solutions exist for 
-Ob < < tt + 9 b with 8 B = sin -1 (2/5^27). The dashed 
line is r — and the dotted line r = \ j\f2. 



The function F{9) acts as the force acting on each oscil- 
lator tending to pin its frequency to the order parameter 
frequency. The range of possible Q for locked oscillators 
is limited by the range of F corresponding to stable sta- 
tionary solutions. 



B. Stability of individual oscillator locked solution 

The stability of the locked solution for an individual 
oscillator mentioned in the previous section is given by 
linearizing Eqs. (21) about the stationary solutions. The 
linearized equations are 



^sin0 



2ar + ^cos<9 ) Sr (58) 



d t 5r = (f3R cos 9) 86 + (1 - 3r 2 )5r 
The eigenvalues are 

A± = 1 - 2r 2 ± VS 



(59) 
(60) 



Then the determinant of the linear matrix derived from 
Eqs. (62) is 



D 



fiR cost 



M 

dr 

-3r 2 



-(l-3r 2 )§+/3i?cos0-g. 

(64) 

The stationary solution r(0) satisfies Eq. (52) and so 
Eq. (64) can be written 



D 



(l-3r 2 ) 



dF 
~d8' 



since 



dF(S) df(6,r(6)) 



dO 



dO 



df 
09 



df dr 
drd§' 



(65) 



(66) 



Thus a zero eigenvalue occurs at and only at station- 
ary points of F(9) or at r = 1/V3- The latter is where 
r(9) has a vertical tangent (cf. Fig. 9) but always occurs 
outside the range of stable solutions, for which we know 
r > 1/V2. 

The only other possibility for an instability is Re A± = 
0, ImA± 7^ 0. This can occur only at r = 1/V2, and if 
S < 0. 

Another result can be derived: if dF/d9 < and 
r > l/y/3 then the determinant D < 0. This implies 
that the eigenvalues are real (since the product of a com- 
plex conjugate pair is always positive), one positive and 
one negative. Thus a negative slope of F(9) implies in- 
stability. Also the Hopf bifurcation can only occur at 
values of 9 where dFjdB > 0. 

To satisfy the imaginary part of the self-consistency 
condition Eq. (50), the range of phases of locked oscillator 
phases must straddle 9 = 0. The oscillator solution here, 
9 = 0, r = 1 is always stable for positive a, B, since 
here A + = -1 + yjl - B(B + 2a). Thus the range of 
possible stable stationary solutions for locked oscillators 
is given by the range of 9 bounded on either side of 9 = 
by the closest stationary bifurcation point or by a Hopf 
bifurcation occurring at r = l/y/2, whichever is closest. 
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C. Properties of the locking force 

We now derive the properties of the locking force, F(&). 
For small B = (3R, Eq. (26) yields a stationary solution 
with r(9) ~ 1 and so 



F(9) ~ B(asin0- cosf 



B^a 2 + 1 sin(0 - 9 a ) 



(67) 
(68) 



with 6 a = tan _1 (l/a). Thus F(9) is a sinusoidal function 
of angle in this limit. The stable solutions are the positive 
slope region for 9 between 6 a — w/ 2 and 9 a + ir/2. 

As B increases, the behavior becomes quite compli- 
cated, and we have not succeeded in proving any results 
about the full range of possible behavior of F(9). Some 
examples are shown in Fig. 10. For B < 2/y/%7 we have 
not uncovered parameters leading to an _F-curve qualita- 
tively different from that in panel (a), i.e. a single max- 
imum and minimum. Note however, that for B > l/v8 
values of r < l/y/2 are encountered, and so the Hopf bi- 
furcation may limit the range of stable solutions moving 
away from 9=0, before the maximum or minimum of 
F is reached. For B > 2/^/27 the physical solutions are 
limited to the range —9b<9<tt + 9b, and the slope of 
r(0) diverges at the end points: dr/dQ — > oo for 9 — ► — 9b 
and dr/d9 — > — oo for 9 — > -n + 6 B . Since 



dF 
d9 



B dv B — — 

— — (a sin 9 — cos 6) H (a cos 9 + sin 9) (69) 

r l d9 r 



we see that when dr/d0 — ► oo 



dF 

~dl 



Fdr 
~d§ 



(70) 



so that dF/d9 — > ±oo with the sign given by the sign of 
F for 9 — > 7r + 9b (where dr/d9 < 0) and opposite to the 
sign of F for 9 — > — 9 B (where dr/dO > 0). This shows, 
for example, that for the parameters such as those in 
panel (c) of Fig. 10 where F(6 — > 7r + 9b) > 0, the slope 
dF/d9 approaches +00 as — ► 7r + 9b- This implies 
that either there is an additional minimum between the 
maximum of F and 9 — > 7r + 0# as in this panel, or the 
maximum disappears and -F(0) becomes monotonically 
increasing in this region, as for the parameters in panel 
(b). On the other hand for a = l,B = 0.5, parameters 
between those of panels (a) and (c), it turns out that 
F(n + 9 B ) < 0, so that dF/d0 — ► -00 as 9 —> n + 9 B , and 
F{9) may decrease monotonically between the maximum 
and 9 — > 11 + 9b- 

In Fig. 11 the regions where the first instability on ei- 
ther side of 9 = is a Hopf bifurcation (Im AqT 7^ when 
Re = 0, with A^ the eigenvalue with larger real value 
for 9 < 0, and A+ for 9 > 0) are plotted as a function of 
a and B. Note that there may be discontinuous jumps 
from Im A 



+ 



to a finite nonzero value of Im A+ : for 
example on the negative 9 side, the minimum in F, giv- 
ing a stationary bifurcation, may disappear by colliding 



with the maximum, and then the first instability jumps 
to the Hopf bifurcation that was previously not the clos- 
est bifurcation to 9 = 0. Figure 11 does not tell the whole 
story about whether the boundary of the locked state is 
a stationary or Hopf bifurcation, because in general only 
the instability on one side of 9 = determines the bound- 
ary, and which side this is depends on the distribution of 
oscillator frequencies via the transverse self-consistency 
condition. This is considered further below. 



D. Self consistency condition 

We have found a range of 9 straddling 9 = giving sta- 
ble locked solutions for individual oscillators. For a fully 
locked solution, all the oscillators must have solutions to 
Eq. (56) in this range. In addition the imaginary part of 
the self consistency condition Eq. (24) must be satisfied. 
This can be written in the form 



g(to) r(tu) sm[9(uj)] du; = 0, 



(71) 



since now there is a unique 9 and r for each oscillator 
frequency. Changing the integration variable in Eq. (71) 
to the angle 9 yields 



g(F(6)) r{9) sin ( 



dF 
~d9 



d9 = 0. 



(72) 



The degree of freedom, for fixed a, B and width of fre- 
quency distribution w, to be determined by this condition 
is the order parameter frequency Q. The scheme of im- 
posing the condition is sketched in Fig. 12. The thick 
solid line, of fixed length w along the F axis, is to be slid 
along the curve F(9) (corresponding to varying 0) until 
the integral Eq. (72) is zero. All quantities, except sin9 
are positive, and so this line must straddle the origin. 
This must be done with all of the solid line lying within 
the stable (dashed) range of F{9). 

For a bounded frequency distribution we can always 
find a fully locked solution for large enough af3. The 
argument is as follows. First, there is always a range of 
stable locked oscillator solutions straddling 6 = 0. If we 
define to = as the center of the distribution g(u>) then 
the order parameter frequency is given by 



-F c -a-B 



(73) 



where F c is F(9) evaluated at the center (with respect to 
the ordinate) of the solid portion of the curve in Fig. 12, 
once Eq. (72) is satisfied. Since it follows from Eqs. (57) 
and (26) that the slope of F(9) at 9 = is 



F'{0) = B(a + B/2), 



(74) 



the range of phase angles of the locked oscillators of order 
w/F'(0) is small for large enough B or aB, and there is 
always a fully locked solution for a bounded frequency 
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FIG. 10: Behavior of the locking force F(8) and the stability eigenvalues X± for (a) a = 1.0, B = 0.2, (b) a = 0.2, B = 0.8, and 
(c) a = 1.0, B = 1.2. The range of angles plotted is 8 between — 7r/2 and 37r/2 for B < 2/^27, panel (a), and between —8b 
and iy + 8b for B > 2/V27, panels (b) and (c). In the eigenvalue plots the solid curves are Re X± and the dashed curves ImA±. 
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_ lmX<t0 / 


Im X> 1 
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lml>= 
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a 



FIG. 11: Regions of Hopf bifurcation (ImA+ 7^ 0) for the 
first instability of the locked oscillator solution either side of 
8 = 0. For 8 < the value of ImA+ is nonzero above the 
dashed curve, and for 8 > it is nonzero above the solid 
curve. 




FIG. 12: Schematic of transverse self consistency condition, 
showing the locking force F(8) acting on each oscillator as a 
function of the phase relative to the order parameter 8. The 
dashed portions of the curve correspond to stable solutions, 
and the dotted to unstable. The thick solid portion denotes 
the distribution of oscillator frequencies that must be placed 
on the stable portion of F and must also satisfy the transverse 
consistency condition. 



distribution in this limit. The center of the band can 
then be evaluated as uj = 0, and so 

= — F(0) — a — B. (75) 

Also, in this limit all the oscillators have essentially the 
same phase, so that R = 1 and /3 = B. This gives the 
complete solution for the fully locked state for very large 
[3. 

It is easiest to understand the limit of the fully locked 
solution by decreasing B at fixed a, w. As B decreases, 



the range of the locking force F decreases. We can con- 
tinue to construct the solution as in Fig. 12, with the 
portion of the F(6) covered by the locked oscillators de- 
termined by the transverse self consistency condition, un- 
til the first value B = B c is reached for which either the 

(s) 

lower end of the locked band would pass below F^ n , 
or the upper end of the locked band would pass beyond 

( s) 

Frn&x- This signals the onset of instability of the corre- 
sponding locked oscillator, either by a stationary or Hopf 
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bifurcation depending on Im A + at the appropriate F u 
or F max [15]. 



This is constructed from Fig. 11, giving the conditions 



0.50 



B 



0.45 



0.35 



0.30 



0.25 




for the instabilities at F r 



is) 



and F^ 



to be Hopf or sta- 



tionary, and the result just determined for which of these 
instabilities limits the range of the fully locked solution. 
Numerical results show that for the top-hat distribution 
the condition /stable < occurs only for a restricted range 
of parameters: large a and B near 2/v / 27, the region to 
the right of the dashed line in Fig. 13. Combining these 
results yields a Hopf bifurcation from the fully locked 
state in the shaded region of Fig. 13. For other oscillator 
distribution shapes we do not have a criterion for the na- 
ture of the instability without a detailed solution of the 
self consistency condition for each width. 

So far the solution has been developed in terms of a, B. 
We now determine the magnitude R of the order param- 
eter from the parallel self-consistency condition, the real 
part of Eq. (23) which can be written in the form 



FIG. 13: Plot showing regions in the B — a plane for which 
the instability from the fully locked state is Hopf (shaded) or 
stationary (unshaded) for a top-hat distribution. The dashed 
line shows the range of a, B for which /stable < for a top- 
hat distribution. (Note the change of range of a, B plotted 
compared with Fig. 11.) 



g{F{6))r{0) cos( 



dF 
d9 



d6 = R, 



(79) 



We could also imagine increasing the width of the dis- 
tribution w from the fully locked solution which occurs 
at small w for fixed a and B, until one end of the grow- 
ing band of locked oscillator solutions reaches F m i„ or 
Fmax- For a top-hat distribution of oscillator frequencies 
the integral 



'stable 



g(F(6))r(6) sm9 



stable 



dF 



d0. 



(76) 



where the integration extends over the whole stable range 
of B, provides an indicator of whether the limit is reached 
at the lower or upper bound of the stable band: if /stable 
is positive, then to satisfy Eq. (72) the integral must be 
reduced by lowering the upper integration bound, and so 
the range of integration extends from the lower stability 
bound. This gives the condition for the maximum w 



w c /2 = F C -F 1 



0) 



(77) 



(s) 

with F±,iL the lowest F for stable solutions and F c the 

mm 

value of F at the center of the band for /stable = 0. On 
the other hand if /stable is negative, then the integral 
Eq. (72) must be increased by raising the lower integra- 
tion bound, and so the range of integration extends to 
the upper stability bound. This gives the condition 



w c /2 = F t 



(s) 



F r . 



(78) 



( s) 

with Flnax the largest F for stable solutions. We can then 
find the range of a, B for which the locked state disap- 
pears by a Hopf or by a stationary bifurcation, Fig. 13. 



where the range of integration is that determined from 
the transverse self-consistency condition, see Fig. 12. 
From R(a, B) we can calculate (3 = B/R at the boundary 
of locked solutions. If the dependence R{B) is smooth 
and monotonic, we can map the results depending on 
{a,B) onto functions of (a, /?). However, discontinuities 
in R(B) might well occur due to the jumps in the sta- 
bility range, for example when the stationary bifurcation 
disappears as described before. This might lead to values 
of (3 for which no prediction, e.g. for w c (a,P), has been 
yielded by the algorithm. We do not have results for such 
cases. 

Two examples of the construction of locked solutions 
for the triangular distribution, defined below in Eqs. (95), 
are shown in Fig. 14: the first for small B = (3R, where a 
solution r(8) exists for all and the second for larger B 
where there are ranges of 8 for which no physical solution 
for the magnitude (i.e. real and positive) exists. The 
value of B for which Eq. (50) is satisfied in each case was 
found by simple bisection applied to the numerical results 
of the integration. In both the cases shown the boundary 

(s) 

of the fully locked region occurs when w m i n = F\^. For 
the smaller value of B this condition corresponds to the 
limit of the existence range of solutions. For the larger 
value of B the limit occurs at the instability of the locked 
oscillator solution, which develops via a Hopf bifurcation 
(Im A+ t^O). Results for the boundary of the fully locked 
state for a top-hat distribution with g(0) = 1 are shown in 
Fig. 19, and for the triangular distribution with g(0) = 1 
in Fig. 22. 
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FIG. 14: Plots of F{&) for a = 2, B = 0.34, w = 1.3 (top 
panel), and a = 2, B = 0.56, u) = 2 (lower panel) correspond- 
ing to the boundary of the fully locked solution: solid curve 
- stable solutions; dashed curve - unstable solutions. The 
dash-dotted lines straddle the range of u = F(6) of locked os- 
cillators, between F m i n and F m i n + w. The values of the order 
parameter are 0.92 (upper panel) and 0.93 (lower panel). 



VII. RESULTS FOR VARIOUS FREQUENCY 
DISTRIBUTIONS 



A. Lorentzian Distribution 

In this section we present detailed results for the case of 
a Lorentzian distribution of frequencies. We concentrate 
on a distribution with g(0) — 1, but also present some 
results for a wider distribution g(0) = |, which shows 
some novel features. We begin this section by calculating 
the linear stability of the unsynchronized state. Since the 
Lorentzian distribution is unbounded, there is no fully 
locked solution. 

For the purposes of analysis we choose, without loss of 
generality, a Lorentzian distribution g(tu) centered about 
zero frequency 



with 



w = MO))" 1 . (80) 



The half width at half height is w. In terms of the shifted 
frequency u) the distribution is 



m = .9(0) 



(Q + 8) 2 + w 2 



with 



8 = a + /3 + 0. 



(81) 



(82) 



The integrals Eqs. (42) are 

2 + w 



h 

h 

h 
h 



2{A + Aw + w 2 +5 2 )' 
8 



4 + Aw + w 2 + S 2 ' 
w 



6 2 ' 



(83a) 

(83b) 

(83c) 
(83d) 



The imaginary part of the self consistency condition 
Eq. (43) reduces to 



2aS + 2w + w 2 



0. 



(84) 



This serves to fix the frequency of the order parameter 
SI at the onset of synchronization via Eq. (82) in terms 
of the parameters of the system a, (3, w. There are two 
solutions for <5 



± \J a 2 - (2w + w 2 ). 



(85) 



For large a or small w the approximate solutions are 

5 ~ (w + T;W 2 )/a giving a locking frequency near the 
center of the band of the free running oscillators, and 

6 ~ 2a giving a locking frequency far in the tails. Note 
that the requirement that 5 is real means that a must be 
sufficiently large |a| > a m in with 



V 7 



2w. 



(86) 



For < \a\ < a m ; n the unsynchronized state is linearly 
stable for all values of f3. 

The critical value of (3 is determined from Eq. (44) and 
evaluates to 



(w 2 + 5 2 ) (4 + 4 w + w 2 + S 2 ) 
(2 w + w 2 ) (a + 5) + 8 (2 - a 8 + 5 2 ) ' 



(87) 



where the expression Eq. (85) for 8 is to be substituted. 
Given a width w of the oscillator distribution, for each 
a > a m i n (w) there are two critical values of /3: (3 C - and 
[3 C+ [corresponding to the minus and plus signs in the 
expression Eq. (85) for 5], such that the unsynchronized 
state is unstable for (3— < (3 < (3+. It is remarkable 
that for very strong coupling, (3 large, the unsynchro- 
nized state remains stable. However, as we have already 
seen and will discuss in more detail, a large amplitude 
synchronized state is also stable in this regime. For 
a ^> a m i n the results for (3 C reduces to /3 C — w(w + 2)a~ 1 
and 4a. In the limit w 1 the former result reproduces 
the result a(3 c — > 2/ng(Q) expected from the reduction to 
the phase equation valid in this limit. 

We have used simulations to confirm the boundaries for 
instability of the unsynchronized state as well as to study 
the behavior subsequent to the instability. The numeri- 
cal results use the cutoff Lorentzian form introduced in 
Eq. (17). Two widths are considered: a narrow one with 
a peak height of <?(0) = 1 and one approximately twice as 
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wide with g(0) = 1/2. In all cases the distribution tails 
are removed above some large frequency. Qualitative dif- 
ferences in the phase diagrams arc observed for the two 
distribution widths. 




P 



FIG. 15: Phase diagram for a Lorentzian distribution of fre- 
quencies with width such that g(0) = 1. Solid and dashed 
lines show analytical results of the linear stability of the un- 
synchronized state. Numerics show the bifurcations are su- 
percritical along the solid portions and subcritical along the 
dashed portion. Dash-dotted lines are saddle-node bifurca- 
tions observed in numerical simulations. States are: U - un- 
synchronized; Si , S$ synchronized with small and large am- 
plitude respectively. 



The phase diagram for the (narrow) Lorentzian dis- 
tribution with g(0) = 1 is shown in Fig. 15. The solid 
and dashed lines are the analytically obtained stability 
boundaries of the unsynchronized solution. The numeri- 
cal simulations show that over the dashed portion of the 
linear instability curve the bifurcation is subcritical, giv- 
ing a jump in the order parameter magnitude R at onset. 
In addition the sweeps yield a number of saddle-node bi- 
furcations identified as discontinuous jumps in R; these 
are denoted by dash-dotted lines in Fig. 15; we do not 
have closed form relations for these boundaries. Thus 
along the dashed or dash-dotted portions of the bound- 
aries, discontinuous jumps in R occur, either between 
the unsynchronized state and a synchronized state, or 
between two synchronized states with different values of 
R. The two synchronized states are labelled Si and S2 
in Fig. 15. When they coexist at the same a and (3 the 
state S2 has the larger value of R, but for both states R 
may go to zero, connecting continuously with the unsyn- 
chronized state U, for some values of a, (3. 

Representative phase diagram slices for the Lorentzian 
case with g(0) = 1, showing the time-averaged order pa- 
rameter magnitudes (R)t as a function of (3 at fixed a, are 
presented in Fig. 16. In agreement with Eq. (86) the un- 
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FIG. 16: Slices of the phase diagram Fig. 15 showing the time 
averaged order parameter magnitude {R)t as a function of /3 
from numerical simulations of the cutoff Lorentzian distribu- 
tion with g(0) = 1. In panel (a) the only solutions are the un- 
synchronized (U) and large R synchronized (S2) states, while 
in panels (b) and (c) the small R synchronized solution is 
also stable. Arrows denote discontinuous jumps that were ob- 
served as simulations followed the various solution branches. 



synchronized state is stable for all (3 provided a < 0.872. 
Simulations find the S2 solution is bistable with U over 
this a range at larger f3 values, as shown in Fig. 16(a). 
For a > a m i n the unsynchronized state is unstable over 
a range in (3. As shown in Fig. 16(b) for a range of a 
near a m i n a subcritical bifurcation occurs at (3 C - as the 
U solution becomes unstable and S2 forms. With in- 
creasing (3 S2 is the only stable solution until a region 
of bistability where both synchronized solutions coexist. 
As shown in the phase diagram there is a small region 
(labelled U + Si + S2) over which all three solutions are 
simultaneously stable. This region of tristability can be 
observed in Fig. 16(b), over the (3 range of hysteric tran- 
sition between S\ and U near (3 = 2.3. With increasing 
a the system passes through a tricritical point where the 
subcritical bifurcation at (3 C - becomes supercritical. The 
tricritical point has a codimension of 2. An example slice 
at a sufficiently large that the bifurcation is supercritical 
is shown in Fig. 16(c), where with increasing (3 regions of 
bistability between the synchronized solutions and then 
U and S2 are observed. 

We now consider a Lorentzian frequency distribution 
of approximately twice the width. Specifically, we take 
g{0) = 1/2 in Eq. (17) or Eq. (80). There are a number 
of changes to the details of the phase diagram that we do 
not discuss in detail. Several of these are straightforward 
consequences of the wider frequency spread: for exam- 
ple the instability of the unsynchronized solution moves 
to larger a and (3. Instead we focus on two particular 
features. 

The first new feature that is evident from the fixed a 
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FIG. 17: Constant a slices from numerical simulations of the 
cutoff wide Lorentzian with g(0) — 1/2. Shown is the time 
averaged order parameter magnitude {R)t for a range of j3 at 
three a values: (a) a = 1.33, (b) a = 1.39, and (c) a = 1.41. 
These solutions are observed in an array of N = 1000 oscilla- 
tors with a Lorentzian frequency distribution with g(0) = 1/2 
and a cutoff frequency of lo c = 16. The dashed box region in 
panel (c) is examined more closely in Fig. 18. 



cuts shown in Fig. 17 is the reconnection of the branches 
of synchronized solutions that occurs as a is decreased. 
In Fig. 17(c), typical of larger values of a, the synchro- 
nized state growing from the /3 C + instability ends at a 
saddle-node bifurcation, with the order parameter jump- 
ing to larger values as (3 is decreased. This is the same 
as the behavior for the narrower distribution, Fig. 16. 
On the other hand for smaller values of a, as in Fig. 
17(a), this state merges continuously with the larger mag- 
nitude state. This change in the topology of the solution 
branches as a increases occurs through the development 
of two additional saddle-nodes, as shown in Fig. 17(b) 
near (3 = 2.1, and the collision of one of these with the 
saddle-node terminating the large-/? upper branch. 

Another apparent difference with the wider frequency 
distribution is the steps displayed by the large amplitude 
synchronized state at larger (3 in Fig. 17. While these 
steps also occur with the narrow frequency distribution, 
their increased amplitude with the wider frequency dis- 
tribution facilitates the presentation of the behavior. As 
an example, the dashed region in Fig. 17(c) is shown 
enlarged in Fig. 18. Panel (a) in this figure plots the 
time averaged order parameter magnitude (R)t, while (b) 
shows that the solution changes correlate exactly with 
changes in the number of oscillators locked to the syn- 
chronized cluster. 

Interestingly, the solution follows different trajectories 
depending on the history of changing (3. For example, 
if (3 is monotonically increased the system settles onto 
the upper trajectory of solution steps (colored red) while 
the lower trajectory of solution steps (colored green) is 
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FIG. 18: Hysteretic behavior between two boundaries of the 
large magnitude {R)t synchronized solution for a Lorentzian 
frequency distribution (color online). Shown are the results 
for N — 1000, a = 1.41, and a cutoff Lorentzian distribution 
with g(0) = 1/2 and uj c = 16. The arrows denote the di- 
rection of the P sweep, and the states are colored red for an 
upward sweep, green for a downward sweep, and black if both 
directions lead to the same state. The + symbols denote the 
starting point for particular sweeps to traverse a state interior 
to the band (green for the start of a downward sweep, red for 
an upward sweep) . Panel (a) is a magnification of the dashed 
box in Fig. 17(c) showing the {R)t over a range of [3. In (b) 
the corresponding size of the synchronized cluster measured 
by the synchronization index is shown. 



observed when (3 is monotonically decreased. These so- 
lutions overlap at smaller (3 (colored black). 

The two solution trajectories are in fact the bound- 
aries of a band of multistable solution branches. This 
structure is demonstrated in Fig. 18(a) by beginning the 
system on each of these two boundaries and reversing 
the direction of (3 variation. Now as (3 changes the so- 
lution moves off its boundary along a continuous line to 
the other boundary where it follows the expected path 
for that sense of (3 changes. Starting solutions are repre- 
sented by the + symbols colored for the direction that (3 
is to be varied. Beginning on the lower boundary, that 
is observed when decreasing /?, and increasing (3 from 
the lower (red) + the oscillators follow the solid (black) 
line in the direction of the upward arrow to the upper 
boundary. Likewise, beginning on the upper boundary, 
that is observed when increasing (3, and decreasing (3 the 
oscillators display a series of states beginning with the 
upper (green) + and follow the solid (black) line to the 
lower boundary in the direction of the downward arrow. 
In this way, the synchronized solution with large time- 
averaged order parameter magnitude turns out to be a 
band of solutions, with striations that can be accessed 
by appropriate changes in (3. The same is true for the 
narrow Lorentzian, however the band is smaller. 

The steps in the magnitude of the order parameter 
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are found to be associated with jumps in the number of 
oscillators locked to the phase of the order parameter. 
Defining a synchronized cluster as the oscillators locked 
to we find the discrete steps in R are due to a small 
group of oscillators leaving the cluster, and then recom- 
bining with the cluster one at a time as (3 increases. In 
Fig. 18(b) the cluster size is measured by the phase syn- 
chronization index, defined by Eq. 15, as the number of 
oscillators with Xn = 1 — e i where e is some small num- 
ber, e << 1, to allow for some phase variation over the 
finite time of measurement. With increasing system size 
any small group of oscillators leaving the cluster can be 
expected to have decreasing influence on the order pa- 
rameter. 

As the number of oscillators tends to infinity the solu- 
tion band becomes a region densely populated by these 
striations as the discrete steps become closer and decrease 
in length. For the simulations in Fig. 18 N = 1000 oscil- 
lators were used. To investigate finite size effects we also 
studied this solution using N = 5000 and N = 10, 000. 
With increasing system size the individual steps on each 
solution boundary occur over a more narrow range in /?, 
becoming shorter in length, and correspondingly more 
densely packed. Thereby, while the discrete steps will 
disappear in the infinite size limit the width of the band 
of synchronized solutions remains nearly the same. In 
this limit synchronized states will move along one of the 
striations to the band boundary appropriate for the di- 
rection of the {3 sweep. 



B. Top-hat Distribution 

For bounded distributions we can calculate the linear 
stability boundaries of both the unsynchronized and the 
fully synchronized states. We first do this for a uniform 
bounded distribution, which we call a top-hat distribu- 
tion. 

A top hat distribution centered on lu = is given by 



w 1 for \oj\ < w/2 
for |u>| > w/2 



The integrals Eqs. (42) for this distribution are 
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£ for 6 < w/2 
for 8 > w/2 ■ 



(89a) 
(89b) 
(89c) 
(89d) 



The self consistency equations (43) and (44) can be solved 
numerically. The equations simplify in the limit of small 
w, and the results here can be displayed in closed form. 



In small w limit there is one solution of Eq. (43) giving 
a locked frequency within the band of (shifted) oscillator 
frequencies, 8 < w/2, for which I\ and 1% may be ne- 
glected in the small w limit. Equation (43) then gives 
the explicit expression for the frequency offset from mid 
band 



c W if 77 

~ — tanh — — 
2 \2a 



(90) 



and Eq. (44) gives the condition of the parameters at the 
onset of synchronization 



2w 



(91) 



Notice that for large a the locking frequency is close to 
the center of the band, and the critical condition reduces 
to 



a/3 c = 2w/tt, 



(92) 



the value expected from the mapping onto the Kuramoto 
model for large a. For small a on the other hand 8 ap- 
proaches w/2 giving us the result that the locking fre- 
quency approaches the upper edge of the band. In this 
case the onset occurs at (3 C ~ 2wa/ir. Even for moderate 
values of a such as a = 1, the locking frequency is far off 
the band center (a fraction 0.92 of the half band width 
for a = 1). 

The second solution to Eq. (43) in the small w limit 
gives an order parameter frequency outside of the band, 
8 > w/2. For small w the frequency is given by 



^coth(f) 
2 \AaJ 



and the onset of synchronization occurs at 



p D ~ 4a. 



(93) 



(94) 



The solution that grows from this instability is a remark- 
able state in that it is synchronized in the sense that there 
is a nonzero value of the order parameter \ip\ ^ 0, but 
there are no oscillators frequency locked to one another or 
to the frequency of the order parameter: a plot of the ac- 
tual frequency distribution as in Figs. 4(a) and 6(a) shows 
a smooth curve with no plateau. Numerical investigation 
of this state shows that instead, the distribution of os- 
cillators p(w, 9) is enhanced for phases near the phase of 
the order parameter: oscillators slow down in this vicin- 
ity (i.e. dO /dt becomes smaller), but do not come to rest 
relative to the order parameter. Now the ordered os- 
cillator frequency plot is continuous, with no plateau of 
locked oscillators. The linear instability boundaries for a 
top-hat distribution of full width w = 1 (giving g(0) = 1) 
are shown as the solid and dashed lines in Fig. 19. The 
approximate solutions Eqs. (90-94) turn out to be quite 
accurate even for w = 1: the approximate curves are 
indistinguishable from the numerical curves in Fig. 19. 
Note that unlike the case of the Lorcntzian distribution, 
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FIG. 19: Phase diagram for a top-hat distribution of frequen- 
cies of width W = 1 such that g(0) = 1. The solid and dashed 
lines are from the linear stability analysis of the unsynchro- 
nized state: the unsynchronized state is unstable for j3 values 
between the two portions of these lines. The numerical re- 
sults show that the solid portion is supercritical, whereas the 
dashed portion is subcritical. Dash-dotted lines are saddle- 
node bifurcations from numerical simulations. The dotted 
line is the stability boundary of the fully locked solution cal- 
culated using the methods of Sec. VI. States are as in Fig. 
15 with, in addition, L denoting the fully locked state. 
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the linear instability persists for a,/3 — ► 0. In this limit 
the frequency of the order parameter is right at the edge 
of the distribution of oscillator frequencies. The discon- 
tinuity of the distribution of oscillator frequencies at the 
band edge seems to be responsible for the persistence of 
the instability to small values of the coupling constants. 



The boundary of the fully-locked solution calculated 
using the analysis of Sec. VI is also shown on Fig. 19 as 
the dotted line. In addition we have performed a care- 
ful numerical scan of the a — (3 plane for N = 1000 or 
10000 oscillators with a uniform distribution of full width 
unity. These results confirm the analytic predictions, and 
again show additional transitions that are inaccessible to 
our analytic calculations. The numerics shows that the 
linear instability from the unsynchronized state is sub- 
critical over the dashed portion of the linear instability 
curves as shown in the figure. Other saddle-node bifur- 
cations are shown as dash-dotted lines. The complete 
phase diagram is quite complicated. Some representa- 
tive numerical sweeps arc shown in Fig. 20. 



FIG. 20: Phase diagram slices observed in simulations with 
a top-hat frequency distribution. The time averaged order 
parameter magnitude {R)t over a range of (3 is shown for 
constant a: (a) a = 0, (b) a = 0.22, (c) a = 0.24, and (d) 
a — 0.60. In these simulations V = 1000 and the width 
w = 1. 



C. Triangular Distribution 

We finally present results for when the oscillators have 
a triangular frequency distribution 



(A/w 2 )(w/2 




\uj\) for M < w/2 
for \uj\ > w/2 



(95) 



This is the case studied in reference [9]. The triangular 
distribution is bounded, and so can have a fully locked 
state as for the top-hat distribution, but does not have 
the discontinuity in g(uj) at the edge of the distribution 
leading to singular behavior as a, (3 — > for that distri- 
bution. The results arc compiled from the linear stability 
analysis of the unsynchronized and fully locked state, as 
well as numerical investigation, usually of 1000 oscilla- 
tors, as shown in Fig. 21. The integrals Eqs. (42) needed 
for the stability analysis of the unsynchronized state can 
again be done in closed form, but the results are too 
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FIG. 21: Solutions observed in simulations of N = 1000 
oscillators having a triangular frequency distribution with 
g(0) — 1. The time-averaged order parameter magnitude 
{R)t is plotted over a range of /3 at fixed a: (a) a = 0.0, 
(b) a = 0.4, and (c) a = 0.9. 



the instability region corresponds to an order parameter 
frequency outside of the band of shifted oscillator fre- 
quencies, so that there is no oscillator with a frequency 
equal to that of the order parameter. This seems to be 
the typical result for a bounded distribution. For the 
unbounded Lorentzian distribution on the other hand, 
there are some oscillators that lock to the order param- 
eter frequency at the large (3 instability. However even 
for this distribution, for most values of a the order pa- 
rameter frequency is far in the tails of the distribution, 
and there are few oscillators at this frequency. Presum- 
ably the effect of these locked oscillators on the critical j3 
is small, so that the dominant physics at this transition 
for bounded and unbounded distributions is not very dif- 
ferent. Near the smallest a where the upper and lower 
stability boundaries meet, the order parameter frequen- 
cies of the two solutions approach the edge of the band 
(from either side). 



VIII. CONCLUSIONS 



cumbersome to list here. 




P 

FIG. 22: Phase diagram for a triangular distribution of fre- 
quencies of full width w = 2 such that g(0) — 1. Symbols and 
lines are as in Fig. 19. 



The results for (3 c (a) for g(Q) = 1 (a width w = 2) 
given by numerically solving Eqs. (43) and (44) are shown 
as the full and dashed lines in Fig. 22. Simulations show 
that this transition is subcritical over the dashed por- 
tions. In addition, as in Fig. 19, saddle-node bifurcations 
identified by jumps in the order parameter magnitude R 
in the simulations are shown by dash-dotted curves. 

As for the top-hat distribution, the large (3 limit of 



In summary, we have analyzed in detail a model for 
the synchronization of nonlinear oscillators where the or- 
dering arises from the reactive coupling between the os- 
cillators, combined with the nonlinear frequency pulling 
of the individual oscillators. Such a model may be a 
more realistic description than previous models for a va- 
riety of physical systems where dissipation plays a rela- 
tively minor role, for example high-Q mechanical oscilla- 
tors and some optical systems. More generally, the model 
may give a more complete description of synchronization 
where the individual frequencies are internally tuned in 
response to a mismatch with other frequencies in the en- 
semble, such as might occur in biological systems. 

We have presented detailed analytic calculations for 
the onset of partial synchronization from the unsynchro- 
nized state, as well as the existence and bounds of the 
fully locked synchronized state at large coupling and non- 
linearity for the cases of bounded frequency distributions. 
The analytical calculations together with numerical sim- 
ulations have been used to construct detailed phase dia- 
grams of the different synchronized states as a function 
of the two parameters of the model, the coefficient of the 
nonlinear frequency pulling a and the coupling constant 
(3, for various frequency distributions. The intersections 
of the various synchronized states lead to rich phase dia- 
grams. 

There are a number of interesting features of these 
phase diagrams. The instability of the unsynchronized 
state occurs over a limited range of (3 at fixed a, so 
that the unsynchronized state regains stability at very 
large values of the coupling strength, although a large 
amplitude synchronized state also occurs here, often ter- 
minated by a saddle-node bifurcation as (3 is decreased. 
This large amplitude synchronized state may also survive 
down to a = 0, i.e. even in the absence of the nonlinear 
frequency pulling. For bounded distributions the state 
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that develops from the linear instability of the unsynchro- 
nized state as (3 is decreased is synchronized (the order 
parameter is nonzero) but there is no frequency locking. 
The phase diagrams also show a wide variety of multi- 
stability, with one or more synchronized states and the 
unsynchronized state coexisting over various parameter 
ranges, leading to hysteresis as the parameters are varied. 
The multistability may be particularly dramatic for the 
large amplitude synchronized states such as displayed in 
Fig. 18 where many synchronized states coexist, leading 
to a band of solutions in the large N limit. 
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